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Abstract 

NASA’s development of new concepts for the Crew Exploration Vehicle Orion presents 
many similar challenges to those worked in the sixties during the Apollo program. However, 
with improved modeling capabilities, new challenges arise. For example, the use of the 
commercial code LS-DYNA 1 , although widely used and accepted in the technical 
community, often involves high-dimensional, time consuming, and computationally intensive 
simulations. Because of the computational cost, these tools are often used to evaluate specific 
conditions and rarely used for statistical analysis. The challenge is to capture what is learned 
from a limited number of LS-DYNA simulations to develop models that allow users to 
conduct interpolation of solutions at a fraction of the computational time. For this problem, 
response surface models are used to predict the system time responses to a water landing as a 
function of capsule speed, direction, attitude, water speed, and water direction. Furthermore, 
these models can also be used to ascertain the adequacy of the design in terms of probability 
measures. 

This paper presents a description of the LS-DYNA model, a brief summary of the response 
surface techniques, the analysis of variance approach used in the sensitivity studies, equations 
used to estimate impact parameters, results showing conditions that might cause injuries, and 
concluding remarks. 


1. Introduction 

During the development of the Apollo command module (CM), many studies were 
conducted to understand the behavior of the CM upon returning to earth. Because of the 
limited computational capabilities at the time, engineers made extensive use of experimental 


1 Trade names and trademarks are used in this report for identification only. Their usage does not 
constitute an official endorsement, either expressed or implied, by the National Aeronautics and Space 
Administration 
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data to complement the relatively simple analyses. These studies were performed to 
understand the landing behavior under a variety of conditions. Interestingly, Apollo was 
initially to include a land-landing architecture. Thus, many of the early studies focused on 
assessing attenuation systems for landing on soil, sand, and clay surfaces [1-3]. However, the 
land-landing capability was abandoned for water splashdown, similar to the Mercury and 
Gemini capsules [4-5]. The exceptions to the water landing studies, were assessments of the 
capsule response in the event of a pad abort [6]. 

The Apollo water landing studies built on the extensive experience of the Mercury and 
Gemini programs. A wide range of experimental parameter studies, using subscale test 
articles, was conducted to evaluate the Apollo capsule performance. A primary interest of 
these studies was to understand the effect of landing attitude and longitudinal or horizontal 
velocity, as well as heatshield flexibility, on capsule response [7-11]. In addition, numerous 
full-scale tests were conducted to assess the capsule performance on more flight-like vehicles, 
e.g., References [12-14], Finally, based on knowledge of probable landing conditions, a 
design envelope was established [15]. The purpose of the present work is to provide a 
probabilistic framework for the Orion Crew Module water landing performance boundary 
evaluation. 

Results from two independent studies are reported here, using data from two sets of LS- 
DYNA simulations, to illustrate the computational approach. One of the simulations used 4 
parameters, whereas the second used 7 parameters to study landing of a rigid capsule on 
water. Although both efforts used the same LS-DYNA model, because the initial conditions 
used in each study are not identical, results cannot be easily combined into a single database. 
Instead results, in the form of time histories, from each study are treated separately and used 
to predict the Brinkley index [16] and to conduct statistical analysis. Fundamental differences 
in the approaches will be pointed out as results are presented. 

For readers not familiar with the Brinkley index, it is a criterion, initially developed for 
emergency escape systems, to assess the probability of injury to a crew member when 
exposed to sudden acceleration pulses, like those observed during landings. Using local 
accelerations at the crew locations in all three directions, the criterion divides the risk of 
injury into three categories: 1) low risk corresponds to a 0.5% probability of injury, 2) 
moderate risk corresponds to a 5% probability of injury and, 3) high risk corresponds to a 
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50% probability of injury. Although there have been efforts to modify the thresholds set 
initially by Brinkley, the calculation of the index is unchanged and is what is used throughout 
the paper. 

The paper presents a computational approach to assess impact conditions (velocities and 
orientation of capsule and impact plane) that could result in injuries to the crew when landing 
a rigid capsule on water. For this problem, a rigid LS-DYNA model is used and described 
first, followed by a detailed discussion of the formulation, the process to create the LS-DYNA 
database, sampling of the parameter space, sensitivity using conventional and analysis of 
variance, development of response surface models, and predictions of injury boundaries using 
the Brinkley index [16]. 


2. Model Description 

The finite element (FE) model of Orion consists of three LS-DYNA [17] parts: the Orion 
capsule, a volume of air, and a volume of water. The capsule is built from Lagrangian rigid 
shell elements using 10307 elements, with an average element size of 3 in. Six nodes in the 
capsule model representing astronaut locations are used for extraction of acceleration results. 
For the water, a total of 1,282,500 cubic elements (4 in. on a side) are used whereas the air 
uses 967,500 elements. Mass and mass moment of inertia properties for the capsule are given 
in Table 1. The volumes of air and water are composed of Eulerian eight-node brick elements 
with dimensions given in Fig. 1. The air is defined with a vacuum material, and the water is 
defined with a null material with a density of 9.3365e-5 lbf s 2 / in 4 and a viscosity of 1 .63e-7. 

The FE model is analyzed using the Arbitrary Lagrangian-Eulerian (ALE) formulation in 
LS-DYNA. The element time step for this model is 22.8 ps, resulting in an analysis time of 
20 hrs per 100 ms of simulation time on a single node of a LINUX-based computer. 
Acceleration results are recorded to the nodout text file. 

Because the impact conditions for the Orion are uncertain, seven critical impact parameters 
have been identified for study, as defined in Table 2, and are assumed to have a uniform 
distribution between the bounds provided in the table. Four of these parameters (horizontal 
velocity, vertical velocity, capsule pitch, and wave slope) are illustrated in Fig. 2. Also noted 
on the capsule are approximate crew locations 1-6. Transfonnations for both the capsule and 
the Euler grid (water and air) are defined in a LS-DYNA keyword file that is separate from 
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the rest of the model. Separate keyword fdes for the Euler grid and the capsule are rotated, 
translated, and linked to the analysis run-stream using the * INC LUDE_TRANS FORM 
option. 

2.1 Mesh Convergence Study 

In this section, the effects of mesh refinement on the FE model are investigated. FE models 
with a rigid capsule at a 38-degree pitch angle and a vertical velocity of 50 ft/s are analyzed 
using 4 different Eulerian models. In all 4 cases, only the thickness of the brick elements in 
the vertical (global X) direction are varied; the in-plane dimensions for the brick elements are 
fixed at 4 in. wide by 4 in. long. Dimensions for the 4 cases are given in Table 3. Also in 
Table 3 are the number of “interface” elements between air and water which are within 8 in. 
of the boundary surface; “far-field” elements are further away from the air/water boundary. 
Results of the mesh convergence study are shown in Fig 3 in tenns of time histories of the 
Brinkley index for cases 1 to 4. Time histories for case 2 are within 0.01% of case 1 (the 
baseline case), but the run time for case 2 is 30% less than for case 1. The maximum Brinkley 
index for cases 3 and 4 are 4.0% and 1.6% higher than the baseline case, respectively. 
Unfortunately, the analysis time for cases 3 and 4 are 4.25 times longer than case 1. Because 
of the computational expense of the more refined analyses (cases 3 and 4), the 4% difference 
in the refined and baseline analyses is considered small, and the solution is considered to be 
converged for case 1 . Hence, results to follow will use the parameters in Case 1 . 

3. Response Surface Description 

Response surface techniques provide a mathematical framework to synthesize models of 
systems for which users have only input/output data. For our purposes, these techniques will 
be used to provide time history predictions of water landings given a limited number of LS- 
DYNA simulations. It is assumed that in order to solve the problem at hand, the amount of 
time to complete all required LS-DYNA runs is prohibited. Therefore an alternate way to 
obtain solutions, i.e., response surface models, is needed. Two techniques have been 
implemented and used; the first one was developed by Sacks [18] is often referred to as 
Kriging and the second technique is the Extended Radial Basis Functions [19-20], Both 
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techniques are briefly discussed to aid the reader in understanding some of the parameter 
settings used in our study. 

3.1 Kriging Response Surface 

In this approach, estimates of the system response are computed using the following 
expression 

Kv) = P + r{v) T R-\Y-p 1) (3.1) 

with 

p = {\ T R l \)~ l \ T R X Y 

m 

R(Vi,Vj) = exp[-0]T (v* -v*) 2 ] 

k=l 

r T (v) = [R(v,v ] i) R(v,v 2 ) ... R(v,v n )] 

y t = [y l y 2 ••• y N ] 

i = 1,2,..., N and j = l,2,...,N 

where vf is k th component of the /' th parameter vector, v e M mxl is a vector representing the 
parameters, N is the population size, y and ft are scalars, R e K /VVA is a correlation matrix for 
the parameters, r(v) e R Wxl is a correlation vector, Y elR^'is a vector with N response 

values, 1 e M^ 1 is a vector whose elements are all ones, 9 is a scalar parameter yet to be 
computed, and m is the number of parameters. Although Eq. (3.1) is written for a system with 
a scalar output, the formulation is applied sequentially for cases with multiple outputs. A 
subtle but critical aspect of this formulation is the fact that the unknown parameter 9 , used to 
define the correlation matrix R , must be computed using maximum likelihood estimation. 
That is, 9 is the value that maximizes the scalar metric 

L(9) = -^Aln(c7 2 ) + ln|i?|), £e[0,oo] (3.2) 

where the standard deviation is computed as 

cj 2 =^{Y-p\) T R\Y-p\) 

Solving this scalar optimization problem produces a response surface fitting error with a 
Gaussian distribution that matches the N response values. However, estimation of the 
parameter 9 requires one to solve a scalar optimization for each output vector; which is 
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conducted here using the genetic optimization algorithm in [21]. For problems containing q 
outputs sampled N times, to generate a Kriging surface requires qN optimizations of the 
metric in Eq.(3.2). Although computationally intensive, the solutions are very accurate. 

3.2 Extended Radial Basis Functions 

The Extended Radial Basis Function (ERBF) method developed by Mullur and Messac in 
[19], approximates the system response as follows 

N N 

y(y) = ^w(\\ v _ v «l) + _ v i) ( 3 - 3 ) 

i=l i = 1 

where y/ is a radial basis function, z i is a scalar parameter solved as part of the least-square 
fitting process, and the function (f> (known as the extended function) takes on different values 
depending upon the value of v-v t , see [19] and [20]. To solve for (/> requires solving a least 
square problem with 3 mN parameters. Since there are N ^values, the total number of 

parameters needed to compute an ERBF model is N(3m+l). In addition, to completely 
define 0, two free parameters must be chosen: 1) the order of a local polynomial (set to 4 in 
this example), and 2) the smoothness parameter y (set to 0.15). Finally, the radial basis 

function y/ is chosen as an exponentially decaying function e~ (v ~ y,) l2, ' e with characteristic 
radius r c set to 0.15. 

In both Kriging and ERBF, scaling of the parameter variables v is extremely important. In 
our implementation the absolute value of the maximum parameter value is used for scaling 
such that when variables are at their maximum value the nonnalized variable equals 1 . 

4. Parameter Sampling 

Parameter sampling is perhaps one of the most important steps when creating response 
surfaces. Techniques such as Latin Hypercube [22], random sampling, Hammersly [23], D- 
optimal[24], Halton-leap [25-26], and even manual setting of parameter variables are just a 
few of the possible means to prescribe a sample parameter. In our study D-optimal, Halton- 
leap, and manual parameter selection have been used at different times for different reasons. 
For example, in the beginning of the study the D-Optimal sampling approach was used 
because it would be followed by an analysis of variance study and this technique is ideally 
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suited for that purpose. On the other hand, for the single variable sensitivity study, being 
conducted in parallel, manual sampling (heuristic approach) made the most sense. 
Nonetheless, for good coverage of the parameter range, Halton-leap deterministic sampling 
approach is also used. 

Later in the paper, results for two independent studies using 4 and 7 independent variables 
are presented. For the 4-variable study, sampling of the parameter variables is conducted 
manually to study specific areas of concern and to facilitate sensitivity studies one parameter 
at a time. In contrast, the 7-parameter study used D-Optimal, Halton-leap, and manual 
sampling. 


5. Analysis of Variance 

Parameter sensitivity in most engineering fields in often associated with derivative 
calculations at specific points in the parameter space. However, for analysis of nonlinear 
systems under uncertainty, sensitivity studies are often conducted using Analysis of Variance 
(ANOVA). In classical ANOVA studies, data is collected from multiple experiments while 
varying all parameters (factors) and also varying one parameter at a time. These results are 
then used to assess the output response variance due to variations of a particular parameter as 
compared to the total output variance when varying all the parameters simultaneously. The 
ratio of these two variance contributions is a direct measure of the parameter importance. 

Sobol in [27] and others [27-31] have studied the problem of global sensitivity analysis 
using variance measures. Sobol developed the "so called" Sobol indices to provide a measure 
of parameter sensitivity. In his formulation the variance is computed using the following 
expressions; 

a = Zt(l) 

1=1 

(5.1) 

i= 1 
N 

D x + M 2 =Zt(l)KA’ z ,) 

(=1 

where // is the mean value, D is the total variance computed using N samples of the function, 
and D, is the single parameter variance due to parameter x. Before using Eq. (5.1) it is 
important to point out a notation change. As before, ia refers to the 7 th sample of the 
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parameter vector but jq is a parameter within v. about which the variance is being evaluated 
and Zi are all other parameters that comprise v. not including x,. To properly evaluate Eq. 

(5.1), one needs at least 2N function (output) evaluations; one where all variables are 
randomly sampled, and a second set where all but jq are re-sampled. With this information 
the 1 st index in Sobol is computed as S x =DJD. Of course, this is only one of the m possible 
factors. Incidentally, Eq. (5.1) is easily extended to study two or more parameter (factor) 
interaction, as described in [27], simply by adjusting the number of parameter that gets re- 
sampled. 

A final note on ANOVA using Sobol's approach is regards convergence of the variance 
estimates. Although Ref. 27 discussed the asymptotic behavior of the variance estimates 
using Eq. (5.1), for cases where a surrogate model is used instead of LS-DYNA, the variance 
estimates are only as good as our surrogate model. Nonetheless, this approach provides an 
excellent way to rank variable importance when only a limited number of LS-DYNA runs are 
possible. 


6. Computational Framework 

To conduct this statistical study, it is preferable to automate the generation of LS-DYNA 
simulations with different parameter values. For this work the computational framework 
described in [32] and [33] is used based on the commercially available MATLAB 1 program. 
MATLAB scripts are used to modify the LS-DYNA keyword input file automatically, to 
update the parameter values from a known distribution, to control execution of LS-DYNA, 
and to read LS-DYNA output files. Afterwards, results are stored within the MATLAB 
environment to conduct ANOVA and optimization studies. Although comparable response 
surface and variance tools are also available in the LS-Opt [17,34], processing of the time 
history data with the existing MATLAB toolbox provides more freedom. 

As mentioned earlier, a companion set of LS-DYNA solutions while varying 4 parameters 
was created manually, outside the aforementioned computational approach, using a heuristic 
sampling approach to prescribe parameter values, to execute LS-DYNA, and to create a 
solution set. These results are also loaded into MATLAB for analysis and are presented later 
in the paper. 
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7. Discussion of Results 


Numerical results are discussed in sections to simplify the presentation and to highlight 
specific issues associated with the various aspects of the problem. First, results for response 
surface (RS) models are presented along with examples of LS-DYNA time history predictions 
using the surrogate model. Following this, results using conventional sensitivity analysis are 
presented. In contrast, results from variance analysis using the approach developed by Sobol 
are discussed using both the 4 and 7 parameter surrogate models. Finally, the RS model is 
used to create contours of constant Brinkley values as a function of impact conditions. 

7.1 Response Surface Modeling Results 

Before RS models are developed, a database must be created with solutions computed 
using LS-DYNA. For this, two independent efforts took place that produced two databases. 
One database contains results for the case when 4 parameters (1- vertical velocity V„ 2- 
horizontal velocity V z , 3- capsule pitch angle 9 V , and 4- wave slope v|/ v ) are varied manually 
(i.e. heuristic sampling), and the second approach were 7 parameters (1- vertical velocity V„ 
2- horizontal velocity V z , 3- capsule pitch angle 0 V , 4- capsule yaw angle 0-, 5- capsule roll 
angle 0 X , 6- Wave slope \\i y , and 7- Wave direction v|/ T ) were varied using a combination of 
sampling techniques which included Halton-leap, D-Optimal, and manual. Each LS-DYNA 
time history output file is read into MATLAB and stored for post-processing. Because the 
initial conditions in the LS-DYNA input key files are set using different logic, the two 
databases cannot be easily combined into one. 

Figure 4 shows results from 231 LS-DYNA runs while varying 7 parameters. To keep the 
amount of data manageable, only 0.2 seconds of simulation time is created and stored, which 
is sufficient to capture the primary impact pulse. At the top plot of Fig. 4a is a summary of 
the time it takes to reach peak acceleration upon impact. This plot is part of the initial 
screening of the database because those cases where the peak time equals 0.2 seconds 
correspond to cases where the capsule is either skimming over the water or cases where the 
horizontal speed is so high that the capsule does not hit the water within 0.2 seconds. Fig. 4b 
shows the resultant maximum acceleration norm in units of gs, and Fig. 4c shows estimates 
for the Brinkley index, computed from predicted local accelerations at the crew location 5. 
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Cases where Brinkley values are greater or equal to 1 correspond to landing conditions with 
probabilities of injuries. 

Once a database is constructed, RS models are created to predict the time responses as a 
function of the input parameters. Figure 5 shows an example of predictions using Kriging 
response surface to predict velocities; solid blue lines correspond to two LS-DYNA solutions 
and the solid green lines correspond to response surface predictions. Admittedly, this plot is 
only a qualitative comparison of the ability of the surrogate model to predict the time 
responses but it shows that the response prediction is within the LS-DYNA bounds. 

A more quantitative comparison illustrating the adequacy of the RS model is shown in Figs. 
6a and 6b where two LS-DYNA solutions are removed from the 4-parameter database and 
time histories are predicted with the resulting RS model. The parameter vector removed for 
each case is shown with the figure. Figure 6a and 5b shows LS-DYNA results in solid blue 
and RS predictions in dashed green. In Figure 6a the error difference between LS-DYNA and 
the RS model is 0.002; whereas, in Fig. 6b, after removing a different parameter set, the error 
difference of 0.03. Although, these are excellent results for these two cases, as is often the 
case with nonlinear problems, there might be regions of the solution space where the 
interpolation is not as good. 

A note regarding the use of Kriging and/or ERBF is in order. Although both approaches 
provide excellent and comparable results, in cases like our 4-parameter case where sampling 
of the parameter space is heuristic, Kriging has problems when the sampling produces ill 
conditioned R(v i ,v j ) in Eq.(3.1). Results with the 4-parameter database are obtained 

exclusively using ERBF, but with the 7-parameters database, both Kriging and ERBF provide 
practically identical results. 

7.2 Conventional Sensitivity Analysis of Brinkley Index 

A standard approach in design is to perform a parameter study, in which one parameter at a 
time is varied while the other parameters are fixed. This approach is sometimes useful to 
show the sensitivity of responses to a given parameter. However, it is highly dependent on the 
selection of the nominal values. As shown in Fig. 7, the Brinkley index is a non-linear 
function of the relative pitch angle (capsule pitch angle minus wave slope). In addition, the 
shape of the Brinkley index curve is different at different combinations of tangential and 
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normal velocity. Clearly, a simple polynomial curve fits will not be able to predict the 
Brinkley index for Orion; more advanced response surface methods must be used. 

7.3 Analysis of Variance 

A more general approach to study the sensitivity of nonlinear system responses to 
parameter changes is the approach developed by Sobol. Normally, Sobol indices are 
computed by direct evaluation of the function values (i.e. peak accelerations, Brinkley values, 
etc) as the parameters are varied. For example, if one needs to evaluate the sensitivity of the 
peak acceleration to variations in pitch angle, LS-DYNA would have to be executed multiple 
times as required by Eq. (5.1). This approach is time consuming with models whose 
simulation times are between 12 to 24 hours. Instead, the RS model is used to provide 
predictions for the Sobol formula. 

To aid in understanding the process, Fig. 8 shows a notional chart where the abscissa and 
ordinate contains Brinkley values for two crew locations 1 and 2. Each dot in the chart 
corresponds to a possible LS-DYNA solution when all the parameters are varied. The outer 
blue circle represents to the total variance, whereas the single parameter contribution is 
depicted using the inner black circle. Finally, the ratio of these two variances provides for a 
direct measure of the parameter importance. 

7.4 Four Parameter Variance Analysis 

To demonstrate the approach for the capsule impacting the water, the 4-parameter response 
surface model is used to estimate the parameter variance using Eq. (5.1). Figure 9 shows the 
fractional contribution from each parameter variation to the total variance of the Brinkley 
index. This fractional contribution is stacked as a percentage of the total variance 
(normalized to 1) observed at each of the six crew locations. Colors are used to distinguish 
among the different parameters, as indicated by the color-bar. Across the top of Fig. 9 are the 
computed mean values for the Brinkley index at each crew location. After reviewing Fig. 9, it 
is clear that the vertical velocity V x is the most important parameter followed by the wave 
slope \|/ x . 
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7.5 Seven Parameter Variance Analysis 

In contrast to the variance results using the 4-parameter surrogate, Fig. 10 shows the results 
using a 7-parameter model. The results are significantly different from the 4-parameter case 
and show that in this case both vertical and horizontal velocities play a significant role along 
with the wave slope. After further review, it is observed that yaw 0- and pitch 0 V show an 
increased contribution to the variance, when compared to the 4 parameter results, and 
significantly change the capsule behavior on impact. 

To better understand these results, it is important to highlight the difference between the 4- 
parameter and 7-parameter databases. Because these were two independent efforts, the LS- 
DYNA initial condition set up is different. The main difference is that in the 4-parameter 
database, the capsule is always placed such that the closest distance from the capsule to the 
water is the same. Hence, during free fall, the capsule always travels the same distance. 
However, in the 7-parameter database this is not the case. Thus results need to be adjusted 
based on impact attitude and speed as opposed to the initial conditions in the LS-DYNA key 
file. 

7.6 Brinkley Contours 

Water impact analysis is aimed at estimating the risk of injury to the crew upon landing. 
For this, the Brinkley index is an industry accepted criteria to assess the risk of injury using 
local acceleration responses at the crew location. If only global outputs are available, the 
method discussed in Appendix A to compute the local accelerations should be used. This 
method is used with the 7-parameter database. 

Brinkley contours are created from the local acceleration for the 4-parameter database and 
the 7-parameter database and compared as a function of the normal and tangential velocity on 
the impact plane. By using the normal and tangential velocities, as opposed to the initial 
velocities prescribed in the LS-DYNA key file, results for the two databases can now be 
compared. Figure 11 shows Brinkley contours at crew location 5 obtained with the 4- 
parameter database shown in Fig. 11a (146 LS-DYNA runs) and the 7-parameter database 
(239 LS-DYNA runs) shown in Fig. lib. Although there are some visible differences in the 
contour curvature and Brinkley values near the extreme velocities, as more LS-DYNA 
solutions are added to the 7-parameter database, results should move closer to those obtained 
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with the 4-parameter data. This is because the 7-parameter RS model requires more LS- 
DYNA solutions to properly capture the behavior of the system when compared to the 4- 
parameter model. Unfortunately the question of how many LS-DYNA runs are needed does 
not have a simple answer. At best, users are encouraged to add solutions in those critical 
regions to see if there are significant changes in the results. 

Brinkley contours provide engineers with a portrait of those conditions that could cause 
injuries. To perform this assessment using a conventional approach, one would need to 
compute enough LS-DYNA solutions to populate the contour. However, this approach would 
only give the user a limited view of the functional relationship for the particular combination 
of variables. In contrast, by using the RS models, a more complete portrait of Brinkley values 
can be generated in which all the parameters are considered. For example, Figure 12 shows 
the Brinkley contours at crew location 5 for pitch angles of -28 degrees Fig. 12a (nominal), 
-18 degrees Fig. 12b, and -38 degree Fig. 12c. Aside from the fundamental differences in the 
curvatures, it should be apparent that at -18 degrees pitch there are many normal/tangential 
velocities that result in Brinkley values of 1 . In contrast, results for the case when the wave 
slope is varied instead between -30 and 30 degrees (with the pitch angle held at -28 degrees) 
resulted in the Brinkley contours shown in Fig. 13. After examining Figs. 12 and 13, it clear 
that changes in the capsule pitch angles are not equivalent to changes in the wave slope angle. 
For completeness Figs. 12 and 13 show the computed impact angles, as defined in Appendix 
B, along with both pitch and the wave slope angles. Note that for a wave slope of 0 and a 
nominal capsule pitch angle of 28 degrees, the relative impact angle (angle between the water 
plane nonnal and capsule impact plane normal) is 28 degrees. On the other hand, a 30 
degrees wave slope corresponds to a 58 degrees impact angle. When studying Fig. 13, 
remember that the independent variables are the initial vertical and horizontal velocities and 
not the tangential and normal velocities. Hence, the normal and tangential velocities are 
derived quantities and therefore certain values are outside the range. 

7.7 Estimation of Critical Impact Conditions 

A more important use of the RS model is to recover critical landing conditions. In Fig. 14 
is an illustration showing Brinkley contours with a fictitious probability distribution function 
superimposed for the normal velocity. If the probability distribution function is known, then it 
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is possible to compute the most probable set of normal velocities that would cause Brinkley to 
exceed 1. This kind of information is precisely what the first and second order reliability 
methods provide. 

To demonstrate the use of RS models to estimate critical landing conditions, the 4- 
parameter RS model is used with a genetic optimization algorithm to compute critical 
conditions. Because parameter distribution functions are not readily available, the optimizer 
is instructed to look for parameter sets that produce Brinkley values of 1. Table 4 shows one 
parameter set of the many sets computed by the optimizer. To verify that this is indeed a 
critical landing condition, LS-DYNA is executed with the parameters shown in the table and 
the time history from LS-DYNA is plotted in Fig. 15 in solid blue versus the prediction from 
the response surface model in dashed green. Although the RS prediction is not exact, it is very 
good, and more importantly it predicted a Brinkley index near 1 . 

8. Concluding Remarks 

The paper discussed an approach that combines solutions from the commercial code LS- 
DYNA with response surface techniques to estimate critical impact conditions of a capsule 
while landing on water. What makes this problem so challenging is the complexity of the 
water and capsule interaction, and the computational time it takes to generate enough 
solutions to understand the physics of the problem. Oftentimes, statistical analysis of 
problems this complex is not undertaken because the computational burden is significant. It is 
shown how to use a limited number of high fidelity LS-DYNA solutions with response 
surface techniques, like Extended Radial Basis Functions and Kriging, to conduct statistical 
analysis and to study conditions under which crew injuries are likely to occur. In addition, the 
approach also allows for sensitivity studies, using an approach developed by Sobol, to screen 
and rank parameters based on their contribution to the output. 

Initial results for the capsule landing on water show that the vertical velocity and the impact 
plane slope are the two most important parameters. Results also showed that although certain 
angle combinations, e.g. capsule pitch and wave slope, produce the same impact angles, the 
system response and Brinkley estimates can be significantly different. 
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11. Appendix A: Estimation of Response for Arbitrary points on Rigid Body 

With the complexity associated with LS-DYNA models, particularly when modeling 
fluids/structure interaction, the amount of data that is created can be overwhelming. For 
storage purposes, it is more efficient to store only a subset of the response data if other data 
can be recovered from it. The ability to recover infonnation from a subset of data is also 
important in those cases where engineers, after completing a set of simulations, realized that 
they should have monitored other points not included in the stored output. Of course, one 
option is to re-run all simulations to store the data for the points left out, but a more efficient 
way, for rigid structures, is to store the least amount of information and then compute all other 
information from the stored data. This approach is, in fact, possible and more efficient 
provided that the initial infonnation stored is judiciously selected to allow for general 3- 
dimensional motion to be recovered. In the following, the reader is reminded of the equations 
to recover the responses at arbitrary points. 

Recall that from rigid body kinematics, the relative velocity of point c with respect to point 
m is given as 


v —v = v , = cox R , 

c m dm dm 


(11.1) 


where cb is the angular velocity of the body and R c/m is the relative position of c with respect 

to m. There are at least two ways that Eq. (11.1) can be used; 1) if the velocity of point c and 
m are known, the angular velocity of the body can be estimated knowing the body’s geometry, 
and/or 2) if the information for an arbitrary point d located at R d/m is needed, given the body 

angular velocity, the relative velocity of d is also given by Eq. (11.1). 

Similarly, if acceleration information is needed, the acceleration of point c with respect to 
m is given as 

A = An + Aim = An + ® X Aim + « X Atm ( 1 1 - 2 ) 


where now the body angular acceleration vector a is introduced. As before, depending upon 
what infonnation is stored with the rigid body simulation, accelerations for points anywhere 
in the body are easily computed using Eq. (1 1.2). 
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For results in the present paper, the angular position, angular velocity, and angular 
accelerations are all available. Hence, Eqs. (11.1) and (11.2) are used with the location of 
new points for estimation of acceleration values. 

Although response extraction points on the body may seem futile, for cases where the 
location of the point at which the response is uncertain, this method provides a simple 
mechanism to examine the variation in the system response to establish uncertain bounds. 
Otherwise, new simulations would need to be executed when evaluating new points. 
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12. Appendix B: Estimation of Impact Velocities and Angles 

In this study, many response quantities are computed directly within LS-DYNA but there 
are others like impact angles and velocities that must be derived from the LS-DYNA outputs. 
In the following, the equations used to compute the impact angles and velocities are provided. 

12.1 Estimation of Impact Angles 

Consider, the capsule shown in Figure B.l. Assume that the capsule motion starts with an 
initial velocity and pitch orientation, falls under a gravity and impacts a water surface that is 
defined by a point R on the surface and its unit normal n p . Also consider a capsule with a 

moving reference frame attach to point m, which is located at the intersection of the capsule's 
center line and an imaginary plane containing the impact point IP. In the following, relative 
vector definitions will have a subscript “a lb” to indicate that it is a vector pointing from point 
b to point a, i.e. making the point b the reference. Vectors without an explicit reference are 
defined with respect to the inertial frame. Finally, if a vector quantity is used without the 
upper arrow, it is referring to the vector of components, for example R alb =[x y 2]’ . 

To estimate the impact angles of a capsule with the water, first the impact point IP on the 
water must be located at the intersection of the capsule and water planes. With the aid of Fig. 
B.l let an impact plane on the capsule be defined by its nonnal n c and a center point m 

located at R m . As mentioned earlier, the water impact plane is defined by its normal n p and a 
point p located at R p . Now the vector difference between point c on the capsule impact 
plane (located at R c ) and R m must be orthogonal to n c , hence; 

n c (R e -R m ) = 0 (12.1) 

Similarly, for an arbitrary point on the water plane IP located at R IP , the vector difference 
between it and R must orthogonal to n and thus satisfies the equation 

n T p (R IP -R p ) = 0 (12.2) 

The plane equations (12.1) and (12.2) are written in terms of vector differences of two 
arbitrary points on each of the respective planes. In addition to these two equations, the 
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impact point, if it exists, is located along the direction B = n c xn , i.e., orthogonal to both 

nonnal vectors. For cases where both planes (capsule and water) are rectangular, B points 
along the intersection line, but, if one of the impact planes is circular (which is our case), the 
intersection line reduces to a single point. 

Because the capsule impact plane normal n c is known or can be estimated from the capsule 
displacements, this nonnal can be used to express any point on the capsule impact plane as 

f = Uq = R c —R m (12.3) 

where the matrix U is an 91 3x2 constructed using the basis vectors from null space of n c , i.e., 

any orthogonal vector set perpendicular to n c ; whereas q is an 91 2x1 vector of arbitrary 

parameters to be computed later. Using the same logic, any point on the impact plane can be 
written as 


w=Gz = R IP - R p (12.4) 

where the matrix G is 91 3x2 constructed from the null space of n and z is an 91 2xl vector of 

arbitrary parameters. In Eqs. (12.3) and (12.4) each of the arbitrary vectors q and z has 
dimensions reduced from 3 to 2. 

Since R IP is the impact point, the following vector summation must hold 

R m +f = R ip (12.5) 

and 


R p +w- Rjp 


(12.6) 


Because the impact point has no projection along B , this fact is represented as 


B T (R m+ f) = 0 


(12.7) 


or equivalently, 

B T (R p +w) = 0 (12.8) 

Combining Eqs. (12.3), (12.4), (12.5), (12.6), and (12.8), yields the following set of 4 
equations and four unknowns 


u ?1M=K + M 

o b t g\\z\ \ -B T R p j 


(12.9) 
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Solving for q and z, and substituting back into (12.5)or (12.6) provides the coordinates of the 
impact point R IP . 

Once the impact point location is known, the impact angles can be easily computed. 
Impact angle#, which is independent of the impact point, is computed using the normal 
vectors as 

# = cos ~\n T c n p ) (12.10) 


Impact angle i// is the angle between a traveling direction R clm and the impact point R IP/m 
(note that both vectors are now defined with respect to m ). Using these definitions the angle 
is defined as 


V 


= COS *( . ^ IP/m ^ dm ) 


R IP J \R 


dm 


( 12 . 11 ) 


As is always the case when defining angles, there is arbitrariness in the way the normals 
are defined and, as such, these definitions must be adjusted to conform to a particular sign 
convention. 
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Fig. B. 1 Capsule and water plane parameter description 


12.2 Estimation of Impact Velocities 

Impact plane location and orientation is provided in terms of a point on the plane and its 
unit normal vector, but to compute impact velocities based on the capsule orientation and 
direction of travel, it is necessary to define the impact plane in terms of a local coordinate 
system. Using the impact plane normal unit vector n a set of orthogonal axes can be defined 

in terms of the null space of n p . If the two bases vectors for the null space are represented as 

G = [gj g 2 ] , one can use the vector 2? = n c x n p , which is normal to both planes, as one of 

the three components of the impact plane coordinate system. If the projection of the vector B 
along gj and g 2 is given as a x - g, ° B and a 2 = g 2 ° B (where ° means dot product), the 
traveling direction is defined as that for which a j is equal to zero. Since our bases vectors are 

not necessarily aligned with any particular axis, one can compute the component along the 
tangential direction (traveling direction) as one of these three vectors; 
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Using w = —n x / as the third component of the coordinate system, the coordinate system is 

now defined by the triad —n l w . The negative sign on the impact plane normal is 

there to accommodate a particular sign convention. 

Finally, using these direction cosines, velocities along the normal and tangential 
components are computed using the dot product of the capsule velocities on impact and the 
impact plane triad. 
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Table 1 Mass and inertial 

1 properties of capsule 

Mass, lbj s 2 / in 

Ixx, lbfS 2 in 

Iyy? lbfS 2 in 

Izz, lbfS 2 in 

49.0536 

216,252 

136,104 

177,684 


Table 2. Parameter definitions 


Parameter 

Number 

Parameter 

Name 

Upper 

Bound 

Nominal 

Lower 

Bound 

1 

Vertical Speed 
V x (in/sec) 

600 

360 

120 

2 

Horizontal 
Speed V z (in/sec) 

0 

-480 

-960 

3 

Pitch Rotation 
0y (deg) 

-18 

-28 

-38 

4 

Yaw Rotation 
0 Z (deg) 

10 

0 

-10 

5 

Roll Rotation 
0x (deg) 

30 

0 

-30 

6 

Wave Slope 
Vx (deg) 

30 

0 

-30 

7 

Wave Roll 
Vy (deg) 

30 

0 

-30 
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Table 3. Solid element dimensions for 4 mesh cases 


Mesh 

Case 

Far-Field 
Air, in. 

Interface 
Air, in. 

Interface 
Water, in. 

Far-Field 
Water, in. 

1 

4 

4 

4 

4 

2 

16 

16 

4 

4 

3 

4 

1 

1 

4 

4 

4 

4 

1 

4 


Table 4. Example of estimated critical landing parameters 


Parameter 

Nominal 

Critical 

Landing 

Condition 

Vertical Speed 
If (in/sec) 

360 

501.9 

Horizontal Speed 
if (in/sec) 

-480 

-47.4 

Pitch Rotation 
0 y (deg) 

-28 

-26 

Yaw Rotation 
0z (deg) 

0 

0 

Roll Rotation 
0x (deg) 

0 

0 

Wave Slope 
Vx (deg) 

0 

-5 

Wave Roll 
Vy (deg) 

0 

0 
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Water/Air Dimensions 


Model 

1 

2 

Width, ft. (Y) 

50.0 

50.0 

Length, ft. (Z) 

50.0 

75.0 

Water Depth, ft. (X) 

19.0 

19.0 

Air Depth, ft. (X) 

14.0 

14.3 


Fig. 1 Model dimensions. 
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Fig. 2 Capsule and water orientation definitions with approximate crew location seating. 
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Fig. 3 Mesh sensitivity results using Brinkley Index. 
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Fig. 4 Summary of LS-DYNA response for 7 parameter database. 
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Fig. 5 Example of response surface predictions of LS-DYNA time. 
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0.4 


Brinkley 

Index 



a) Parameters removed P=[360 in/s, -480 in/s, -28 deg,0 deg] 


Index 



b) Parameters removed P=[540 in/s, -480 in/s, -28 deg,0 deg] 


Fig. 6 Verification of response surface and LS-DYNA time history predictions at 
crew location 4 after removing solution from the response surface, fitting 
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Fig. 7 Single variable sensitivity analysis. 
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Fig. 8 Example of analysis of variance. 
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Brinkley Index Mean Values 



Fig. 9 Variance analysis of Brinkley index using the 4 parameter RS model. 
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Fig. 10 Variance of Brinkley index using the 7 parameter ERBF model. 
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Fig. 1 



a) 4-Parameter database and 146 LS-DYNA runs 



b) 7-Parameter database and 239 LS-DYNA runs 
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1 Comparison of contour plots for the Brinkley index using the 4-parameter 
and a 7-parameter database. 
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Fig. 12 Contour plots for Brinkley Index as a function of the impact tangential 
and normal velocity at the crew location 5. 
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Impact angle 9=58 deg 
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Fig. 13 Contour plots for Brinkley index as a function of the impact tangential and 
nonnal velocity at the crew location 5. 
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Crew Location= 1 



Fig. 14 Illustration of process to estimate most likely critical condition. 


Crew Location 4 Predictions 



Fig. 15 Verification of predicted critical landing condition using the 
response surface model. 
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